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We analyze the stability of dust layer in protoplanetary disk to understand 
CLc the effect of the relative motion between gas and dust. The previous analyses 

Q ■ not including the effect of relative motion between gas and dust show that the 

shear-induced turbulence may prevent the dust grains from settling sufficiently to 

■ be gravitationally unstable. We determine the growth rate of Kelvin-Helmholtz 
instability in wide range of parameter space, and propose a possible path toward 

^ ■ the planetesimal formation through the gravitational instability. We expect the 

■ density of dust layer becomes Pd/Pg ~ 100 if the dust grains can grow up to 10m. 

Subject headings: planetesimal formation, turbulence, shear instability, gravita- 
tional instability 



1. Introduction 

The recent discoveries of extrasolar planets are indicating that the formation of Jupiter- 
mass planets is a common process. According to the current standard model of planet 
formation, the growth of planet-sized body occurs in two sequential phases. In the first phase, 
micron-sized dust grains grow up to the kilometer-sized bodies. The kilometer-sized bodies 
called planetesimals are the building blocks of terrestrial planets, cores of giant planets. 



- 2 - 



comets, and asteroids. In the second stage, the planetesimals continue to grow through 
inelastic coUisions. 

The initial stage of coagulation of dust grains is supposed to continue up to about cen- 
timeter size. However the growth of dust grains, from cm to km size, is not well understood. 
The agglomerative growth from submicron-sized to kilometer-sized bodies has many prob- 
lems. The inward orbital drift associated with gas drag that would carry such mctcr-sized 
bodies from 1 AU into the central star is very rapid, and its timescale is only 10^ years 
(Adachi, Hayashi, & Nakazawa 1976; Weidenschilling 1977). The direct assemblage of dust 
grains in laminar flow may be impossible. 

One of the scenarios for this stage is due to gravitational instability (GI). If particles 
settle into a layer having sufficiently high density and low velocity dispersion, high den- 
sity region spontaneously collapses under their self-gravity forming km-sized planetesimals 
(Safronov 1969; Goldreich & Ward 1973; Coradini, Magni, & Federico 1981; Sekiya 1983; 
Yamoto & Sekiya 2004). 

There is a critical issue in this GI scenario. As dust grains settle toward the midplane, 
the rotational velocity around the midplane increases because of the reduced effect of the gas 
pressure compared to the centrifugal force and the solar gravity. The rotational velocity is a 
function of distance from the midplane, and the shear may induce Kelvin-Helmholtz insta- 
bility (KHI). The shghtest amount of turbulence in the nebular gas due to KHI prevents dust 
grains from settling. Many authors have investigated this issue (Weidenschilling 1980; Cuzzi 
et al. 1993; Champney et al. 1995; Weidenschilling 1995; Sekiya 1998; Dobrovolskis et al. 
1999). In the case of the minimum mass solar nebula, they concluded that the gravitational 
fragmentation of the dust layer is impossible if the turbulence is developed. Sekiya (1998) 
showed that large values of the total dust-to-gas mass ratio may provide a density cusp at 
the midplane. Youdin & Shu (2002) and Youdin & Chiang (2004) argued this cusp as a 
triggering of GI. 

In order to understand the effect of instability in detail, linear stability analyses were 
performed. Sekiya & Ishitsu (2000) have performed the linear analysis of the shear insta- 
bility for the constant-Richardson-number dust density distribution. Their result confirms 
standard expectations by showing that the midplane shear layer becomes unstable when the 
minimum value of the Richardson number decreases below 0.226. Sekiya & Ishitsu (2001) 
investigated the instability using the hybrid dust density distribution in which the dust layer 
has a constant density region and transition region. For this distribution, the growth rate of 
the shear instability is much larger than the Kepler angular frequency when the dust density 
at the midplane is larger than the gas density. Ishitsu & Sekiya (2002) investigated the shear 
instability of the hybrid dust density distribution including the Coriolis force but neglecting 
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the tidal force. Their results showed that the Coriohs force has little effect on the growth 
rate of the shear instability. Ishitsu & Sekiya (2003) investigated the shear instability of 
the hybrid dust density distribution including both the Coriolis force and the tidal force. 
They showed that although the tidal force tends to stabilize the shear instability, the shear 
instability occurs before the dust density reaches the critical density of the gravitational 
instability. Throughout those studies they use a single-fluid approximation. Garaud & Lin 
(2004) performed the hnear stability analysis based on two-fluid formahsm with strong cou- 
pling approximation (a < 1cm) and destabilizing effect of radiative cooling. They confirmed 
that the shear instability occurs prior to gravitational instability. 

In this paper, we remove the strong coupling approximation from the previous analyses, 
that is, we take into account the effect of the relative motion, neglecting the Coriolis force 
and the tidal force. In Section 2 we investigate the effect of the relative motion between gas 
and dust in the two uniform fiuids in relative motion separated by a horizontal boundary. We 
derive the condition of stabilization by the relative motion. In Section 4, the basic equations 
for the linear analysis in a realistic protoplanetary disk are derived. In Section 5, results are 
given. In Section 6, we discuss a possible path toward the planetesimal formation by the 
gravitational instability. Section 7 is for conclusion. 
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2. Kelvin-Helmholtz instabilities in the case of two uniform dusty fluids in 

relative horizontal motion 

The linear Kelvin-Helmholtz instability of shear layers is well-known for flows without 
dust grains (Chandrasekhar 1961). We consider the effect of dust in the Kelvin-Helmholtz 
instability. First we examine the dispersion relation for the linear perturbation on Heaviside 
step function velocity profile in this section. Since this velocity profile is very simple, one 
can obtain the analytic dispersion relation, which help us to understand the character of 
this instability. We analyze this dispersion relation to evaluate the stabilizing effect due 
to the inertia of dust. In the following sections we investigate the linear stability of the 
Kelvin-Helmholtz instability in a more realistic protoplanetary disk. 

We suppose that the streaming takes place in the y-direction with the velocity U{z), 
and gas and dust have the same velocity in the unperturbed state: 



U{z) 



[/_ forz<0 
C/+ forz>0 



The unperturbed densities of the gas and the dust is uniform everywhere. The gas is in- 
compressible and inviscid fluid. The velocity dispersion of dust is negligible. For simplicity 
we do not take into account the influence of gravity in the vertical direction and the size 
distribution of dust grains. 

The equations governing this system are 

V-Vg = 0, (2) 

^ + V(pdVd) = 0, (3) 



+ (vg • V)vg^ = -VP, + ^PgPd(vd - vg), 
Pd + (vd • V)vd^ = -Apgpd(vd - Vg), 



(4) 



(5) 



where A ,Pg,Pd, Pg, Vd, and Vg are the friction coefficient, the gas pressure, the dust density, 
the gas density, the dust velocity, and the gas velocity, respectively. We assume the form of 
perturbation quantities as 

q{z)ex.p{-i{u;t-ky)). (6) 

We only consider the perturbed motion in the yz-plane since these modes correspond to 
the most unstable modes of the Kelvin-Helmholtz instabihty. To obtain the equations 
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for the perturbation, we set the quantities U = {0,U{z),0),'Vgi — (t'gxi, 'I'gyi, 'i'gzi),Vdi = 
(^^dxi, -^dyi, -^dzOiVg = U + Vgi,Vd = U + Vdi, and P = Pq{z) + Pi{x, y, z). 

The perturbation equations are 

ikv^,{z) + '^-^{z)^Q, (7) 

-i{uj - kU±)pdi{z) + ikpdoVdyi{z) + Pdo^^(^) = 0, (8) 
-ApgoPdo{vd^i{z) - Vg^i{z)) - ipgoVg^i{z){u; - kU±) = 0, (9) 

-ApgopM{vdyi{z) - Vgyi{z)) - ipgoVgyi{z){uj - kU±) + ikPi{z) = 0, (10) 

dP 

-ApgoPdo{vdzi{z) - Vg^i{z)) - ipgpVg^i{z){uj - kU±) + -^{z) = 0, (11) 

-ApdoPgo(Vi(^) - ^dxi(^)) - ipdoVd^i{z){uJ - kU±) = 0, (12) 
-Apaopgoivgyi{z) - ^;dyi(^)) - iPdoVdyi{z){u; - kU±) = 0, (13) 
-Apdopgo{vg^i{z) - Vdzi(^)) - ipdoVdzi{z){uJ - kU±) = 0. (14) 

We use the x-component of these equations (9) and (12), we have 

uj - kU± = -iA{pgo + pdo)- (15) 

This stable mode corresponds to a simple decay of the relative motion of the dust and the 
gas owing to the friction. 

Eliminating the Vgyi{z) , Vdyi{z) , v^ziiz) , Pi{z) from the equations (7), (10), (11), (13), and 
(14) , we have the equation of Vgzi{z): 

^{z) = ev^,{z). (16) 

To obtain the condition that the perturbed quantities that do not diverge at 2; = ±00, we 
have the solution of the equation (16): 

/ \ _ j C+exp{—kz) for 2; > , . 

Vgzi(^j = | C_exp(A;z) for ^ < ' ^ ^ 

where C+,C_ arc arbitrary constants. 

We need two boundary conditions. First we consider the boundary condition for velocity. 
The z-component of the velocity vector Vgz at the interface is related to the position of the 
boundary Zh{t, x, y): 

dzy^ dz\) dzyy 
.g. = ^ + .g.^ + .g.^. (18) 
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The perturbation equation is given by 



Vgzi(x, y, ±0, t) = U±^{x, y, t) + ^{x, y, t). (19) 



We assume 
and we have 



Zb{x,y,t) = Zhexp{-i{ut - ky)), (20) 

^g^l(^O) 

= ^7777 N ■ (21) 

i[kU± — CO) 

We obtain the boundary condition that Vgzi (z) / {i{kU {z) —a;)) should be continuous at 2; = 0. 
This condition leads to 

^- (22) 



u — kU+ ijj — kU- 

However we do not need this boundary condition for dust (see Appendix A for the reason). 
Second we must assume that the pressure is continuous at {x, y, z\){x, y)). P\{z) is given 



by 



and we obtain 



A;^ Apgo — i[ui — kU±) az 



^ ' Apgo - - kU+) ^ ' Apg^o - i{uj - kU_) ^ ' 



With these conditions we can find the equation of lo by eliminating C. 



± • 



4 .7^d + 2 . 7ed + l 2 iU^k\nd-2) U^eiTZd + l) U^k^ 
r + I a; ^ uj -p -uj 

^drag ^drag 4idrag ^^dv&g 



where 7^d = Pdo/PgO, ^drag = iMPgo, and U± = ±f . 

In the limit of either TZ^ 0, tdrag — >^ 00, or tdrag — >^ 0, we have 

kU , , 

a; = —I, (26) 

which is the simple unstable mode of the original Kelvin-Helmholtz instability. If there is 
no dust (T^d 0), obviously the dispersion relation reduces to that of the original Kelvin- 
Helmholtz instability. If the friction is very small, the gas can move freely without the 
effect of dust, thus it reduces to the original Kelvin-Helmholtz instabilities. Furthermore, 
if the friction is very strong, the dust cannot move independently, that is, gas and dust 
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behave as one strongly coupled fluid, and dispersion relation reduces again to the original 
Kelvin- Helmholtz instability. 

The equation (25) depends on kU. We introduce u; — iCjkU/2 and tdrag — tdrag/ (kU), 
and we obtain the equation that does not depend on kU : 



,4 , 2(7ea + 2) ._3 , A{na + 1) ^_, , 2(7^d - 2) 4(7^d + l)+tLg . 
uj H J u H u H J u ^ U. [27) 

''drag f-drag ''drag f-drag 



a;(fdrag)^d) determine the degree of instability. In the limit of ^drag 0, oo, we have 
a;(tdragj "^d) — 1- This equation has always a real solution and we consider only a real 
solution (see Appendix B for reason). 

In Figure 1, we plot the contour of a)(tdrag, "^d)- The u approaches to 1 in TZ^ <^ 1 or 
^drag — ±oo. It corresponds to the region where the friction is not effective for the reduction 
of growth rate. In the intermediate region, 71^ > 1 and <^ ^drag <^ the growth rate is 
less than that of the original Kelvin-Helmholtz Instability. 

We consider the condition of stabilization. The friction should not be strong or weak to 
stabilize Kelvin-Hclmholtz instabilities. Wc express this condition as min(l/Apg, 1/Apd) < 
1/kU < max(l/Apg, 1/Apd)- We need also that the amount of dust is more than that of gas, 
> 1. 

The growth rate is reduced less than /i {uj < /i) ii the condition is satisfied, which is 

Creducel < ^^^drag < C'reduce2 (28) 

where 



{n^ -2)ii + {na + 2) - •^7^dV'(l + 1^'? - 47^d (i - 1^") - 4(i - i^'f 



Creducel — 14 



(29) 



_ (7^d - 2) /i + (7^d + 2) + ^7^dV(l + /^')' - 47^d (1 - /i^) - 4(1 - /i2)2 

C'reduce2 — ^ ■ (30) 



/i' 

In the case of TZ^ 3> 1, we have the approximated condition: 



< A;C/Wg < T^^d. (31) 



In addition, in the case of <^ 1, we have the simple condition: 



2 

- < /ct/idrag < 2/x7ed. (32) 
A* 
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Figure 1 shows that the growth rate in the region satisfying this condition is reduced. 

We study how the friction reduces the growth rate of the Kelvin-Helmholtz instabihties. 
We seek the minimum of a)(idrag) ^d) in — const.. In the case of TZ^ ^ 1, one can obtain 

N 2 kU . 

^^min(7^d) ^ —7=—l, (33) 

Wg - 2V^^- (34) 

As shown in Figure 2, this expression is in good agreement with the numerical solution for 
7^d > 10. 
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Fig. 1. — The contour of growth rate a; as a function of Pd/Pg and ^drag- 
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Fig. 2. — The maximum stabilized growth rate. The sohd hnes denotes the numerical 
solution, and the dashed lined denotes approximate solution. 
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3. The case of dusty layer between two gaseous layers 

In this section we consider the effect of the thiclcness of the dust layer. The dusty gas 
layer whose thickness 2d is sandwiched between two gas layers without dust. For simplicity 
we assume the unperturbed state is symmetric with respect to the equatorial plane {U {z) = 

= p N < , (35) 
^ ^ \ for 1^1 >(i ' ^ ^ 

. N / Pdo for \z\ < d 
P'^^'^ = [ for |.| > d ' ^''^ 

Pgo(^) = Pgo- (37) 

The perturbation equations are the same in Section 2. Thus we have the perturbed 
velocity that does not diverge at 2; = oo: 

, , r C^e^' + C^e-^' for < z < d 
"-^^^^ = 1 Cae-^^ ford<. ' ^^'^ 

where Ci,C2, and C3 are constant. 

From (21) and (23), we have two boundary conditions at 2; = ^d- 

(2kdr r ^ {kU -uj){kU-uj- iAjpap + pgp)) 

^'-^'^ kU-u;- iAp^ ""^^ - ° ^^^^ 

e^'^'^uCi + CUC2 - (cu - kU)C3 = 0. (40) 

Because of the symmetry of the unperturbed state, there are two types of solutions: even 
solutions where Pi{z) = Pi{—z), Vgzi{z) = —Vgzi{—z) and odd solutions where Pi{z) = 
—Pi{—z), Vgzi{z) = Vgzi{~z). The boundary conditions at z = are written as 

Ci + SC2 = 0, (41) 

where S — ±1. The solutions in the case of 5" = 1 correspond to the even solutions, the 
solutions in the case of 5" = — 1 correspond to the odd solutions 

By the condition that equation (39), equation (40), and equation (41) have a nontrivial 
solution, we obtain the dispersion relation: 

^3 + (7^de-2%/,g5 + (2 + Ua) fiig + 2i {Se-''^' + 2)) /.^ 
+ 2i {e-'''S + 1) (2(1 + 7^d)^ii, + {e-'''S + l) ((1 + 7^d)^iig + i) = 0(42) 



- 12 - 



where = uj/{kUi/2), fdrag = tdiagkU, TZ^ = Pdo/Pgo- 

If the friction is very weak ( t^rag ^ 1 ) ; the growth rate Im(a;) is given by 

HI , 

Im(u;) = ^Vl-e-4fc<^. (43) 

In this case, the odd mode and even mode have the same growth rate. 
If the friction is very strong, the growth rate Im(a;) is given by 



V(i + ^d)(i 



e 



kU \^ ,^ rr-jT — (odd mode) 

Thus odd mode is more unstable than even mode. In the limit of 71^ 0, the growth rate 
for these modes approaches to the growth rate for the case of weak friction. 

In the long wavelength limit {kd — > 0), one can show 

2 



-.{kdy^'^ — 2i (even mode) 



Vi + n~d ' ' , (45) 

Cl {kd) V2 - {kd) V2 (odd mode) 



where c\ and ci are given by 



4(l + 7^d+^Lg) 

''drag 



- - \- (46) 



ciC2 = ~2 • (47) 

"I" ''drag 



Prom the equation (46), we give the inequality: 

.l>.^^^A..-%^>.>' (48) 

This inequality means that the odd mode is always more unstable than the even mode in 
the case of the long wavelength limit. We can obtain the analytic solutions of C\ and Ci'- 



+ ''drag + ''drag 



= 2W^ 4- ^-2 ''"^^ 50 

^ ''drag ^ ^ ''drag 
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In the short wavelength hmit {kd — > oo), /i is constant for k. There is no difference of 
the growth rate between the even mode and the odd mode. 

The growth rate is proportional to k {kd — > oo) or A;^/^ (kd — > 0). In the short wavelength 
limit the growth rate diverges owing to the discontinuity of the unperturbed velocity. In the 
Section 5, we show that the growth rate does not diverge in the short wavelength limit if the 
unperturbed velocity is continuous. 

In Figure 3 we show the dispersion relation in the case of idrag — 0.1C//d, TZ^ — 10. In 
the case oi kd < 1, there is the difference between the even mode and the odd mode, and 
the growth rate is proportional to A;^/^. In the case oi kd > 1, the even mode and the odd 
mode have similar growth rate, and the growth rate is proportional to k. 
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Fig. 3. — The growth rate as functions of wavenumber k for odd mode (sohd curve) and 
even mode (dotted curve) in the case of tdrag = 0-lU/d, TZa = 10. 
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4. Basic equations 

In this section we analyze the stabihty of a structured dust layer in the protoplanetary 
disk, which is more realistic than the step function velocity profile studied in the previous 
section. 

We do not assume that the dust aggregates are small enough to couple firmly: the 
dust aggregates can move relatively to the gas. For simplicity, we use the local Cartesian 
coordinates {x, y, z) and neglect the curvature of the orbit of the perturbed state, and do not 
consider the gravity of central star and the effect of Cohoris force. The equations governing 
this system are 

V • Vg = 0, (51) 



dt 

dv 



dt 



+ V(pdVd) = 0, (52) 
+ (vg • V)vg^ = -VP, + Apgpd(vd - vg), (53) 

Pd + (vd • V)vd^ = -/lpgpd(vd - Vg). (54) 

We assume the gas density pgo is uniform {dpgo/dz = 0) in the dust layer since we 
suppose a thin dust layer. We adopt the model of dust density Pdo{z) which has the form as 



Pdo for \z\ < Zd- 2hd 



Pdo{z) 



PdO 

2 



. z- Zd + hd 
1 — smTT 



for Zd — 2hd < \z\ < Zd , (55) 



2/id 

for Zd < z 



where Zd the half thickness of the dust layer, hd the half thickness of the transition layer, z 
the distance from the midplane, and pdo the dust density at the midplane (Sekiya & Ishitsu 
2001). Zd is given by 

Zd = h hd, (56) 

2pdo 

where Ed — pd{z)dz is the surface density of dust. 

We use here a frame of reference moving with gas dominant region (pdo(^) = 0). The 
stationary rotational velocity of gas is given by Nakagawa ct al. (1986): 



(57) 



-16- 



where TZd{z) — Pdo{z)/pgo, T = pgA/D,-k_, i>k the circular Kepler velocity, Qk Kepler frequency, 
A the friction coefficient. 77 is a non-dimensional parameter which represents the effect of 
the radial pressure gradient: 

2^ai^' ^^^^ 

where P is the gas pressure, Cg the isothermal sound velocity, and r the distance from the 
rotational axis. We assume that the friction coefficients for particles whose sizes are smaller 
than the mean free path of the gas are determined by the Epstein's drag law, and those 
larger than mean free path are determined by the Stokes law: 

for a < -Amfp 



where Amfp is mean free path of the gas, a the radius of dust aggregate, pmat the internal 
mass density of dust aggregates, and Ct the thermal velocity of gas. We set pmat — lgcm~^. 
For the sake of simplicity, we neglect the radial velocity of gas, and assume that the dust 
and the gas have the same velocity in the unperturbed state. The equations are similar to 
those of Section 2, except that the unperturbed velocity and the unperturbed dust density 
is not uniform. We restrict ourselves to the case where = since this gives most unstable 
mode. We have the perturbed equation: 

(Pi' dv 

= Ci{U{z),pM{z),pg,uj,k,A)-^ + C2{U{z),pdo{z),pg,ij,k,A)vg^i. (60) 

Ci and C2 are 

Ci = \' , (61) 

1 + l^d — itdragUJo 
^ _ k{i + tdragl^o) ^ 

UJo{l + TZd- ^^drag<^o) 

(« + ^drag<^o) (« + ^drag<^o) (« + ^drag<^o) / 

where a prime means the derivative with respect to z, 

cuq — u! — kU{z), (63) 

id.ag = (64) 
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In the limit of idrag — this equation approaches to equation (22) of Sekiya & Ishitsu (2000) 
without the effect of gravity. The perturbed quantities are given by 

p ^ Pg0^o(-l " 7^d + itdrag^o) dVg^i ^ ^^'(^)PgO " (i+Wga;o)^) ^ 

^ A;2(i + tdragi^o) dz k 

_ -i^gz/(^) + Wg(fe^gzl(^)^'(^) -UJQVgzl'iz)) , . 

M^ + Wgc^o)^ ' ^ ^ 

v.. = (68) 

In region outside the dust layer z > z^, equation (60) becomes 



= k\,i. (69) 



dz^ 

We find the solution of the perturbed quantity that does not diverge at z — > oo: 

Vgzi (z) — C+ exp(—kz) for z > z^ , (70) 

where C+ is a constant. 

In the case of ^ < 2;d, we cannot obtain the analytic solution. Thus we solve equation 
(60) numerically. 

We consider the boundary conditions. First Vgzi{z)/{u; — kU{z)) should be continuous 

aX, z — z^: 

^gzl(^d + 0) ^ t^gzl(^d - 0) 

uj-kU{zA + Q) u-kU{zd-0)' ^ ' 

Since the unperturbed velocity is continuous, this condition means that Vg^i is continuous at 
z = z^: 

Vgzi {zd + 0) = Vg^x (^d - 0) (72) 

Second the gas pressure should be continuous. With equation (65) , (70), and (72), we obtain 
the boundary condition at z = ^d- 

A;^gzi(^d-0) + ^(^d-0) = 0. (73) 

There are two types of solutions: even solutions where Pi{z) — Pi{—z), Vg^i{z) — 
~'^szi{~z) and odd solutions where Pi{z) — —Pi{—z), Vg^i{z) — Vg2\{—z). We have the 



boundary condition at z—0: 



Vgzi = (even mode) 

^ = (odd mode) ' ^^'^ 
dz 

In the next section we solve equation (60) with the boundary condition (73) and (74) 
numerically. 
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5. Numerical Result 

We adopt the model of the minimum solar nebula at r = lAU (Hayashi 1981): 

7.l(^)-'/'g/cm2 for (0.3AU < r < 2.7AU) 
30 (^)"^^^g/cm2 for (2.7AU < r <36AU) ' 

rj = 1.81 X 10-3 (— j , (76) 

/ T \ —11/4 

Pgo = 1.4 X 10-^ (^) , (77) 
Amfp = 1.4 (j^) cm, (78) 

/ r \ —1/4 

Cs = 1.1 X 10^ j cm-s-^ (79) 



As shown in Figure 4, the hybrid dust density distribution has a constant density region 
and a transition region whose width is 2/id. In the limit of 0, the dust density 

distribution becomes the step function discussed in the Section 3. In Figure 5, we can see 
that the eigenfunction of v^zi approaches to the cigcnfunction in the case of the step function. 
We restrict ourselves to the case of = O.bzd. in the following analysis. The unperturbed 
state is characterized by two parameters: the size of dust a and the ratio of the dust density 
to the gas density Pd/Pg on the midplane. Figure 6 shows the distribution of the unperturbed 
density and unperturbed velocity. 

Figure 7 shows the growth rate as a function of wave number in the case of a = 0.01cm 
and Pd/Pg = 0.1. In the limit of /c — > 0, the growth rate approaches to 0. The growth rate 
has a peak at a finite wavenumber. The modes with short wavelength {k > kc) are stable. 
The odd modes are more unstable than the even modes. Thus we restrict ourselves to the 
odd mode. 

In the hmit of A; — > 0, (i.e., A = 27r/k ^ hd — Zd/2), the width of transition region is 
much smaller than the wavelength, thus we expect that the velocity profile approximately 
corresponds to the step function discussed in Section 3. As shown in Figure 8, the growth 
rate is approximately proportional to k^^^, and approaches to as A; —> 0. 

In the limit of A; — > oo, the wavelength is much shorter than the width of transition 
region. If the unperturbed velocity is continuous, the velocity shear is negligible, hence the 
flow is stable in this limit. 
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Fig. 4. — The unperturbed density and the unperturbed velocity in the case of Pd/Pg = 1, 
a = 1cm. The sohd hne corresponds to step function density profile discussed in Section 
3. The dashed line {h^ = OAz^) and the dotted line (/la = 0.5Zfi) correspond to continuous 
density profile. 




Fig. 5. — Same as Fig. 4, but for the eigenfunction of Vgzi- 
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In Figure 9, the amplitude |-Pi(2;)| and the phase ^{Pi{z)) are drawn in the case of 
a = 1cm and Pd/Pg = 1, where $(/) = tan"-^ Im[/]/Re[/]. |-Pi(^)| have a maximum value 
at z/hd ~ 1.2, and decreases to as z/h^ oo. $(Pi(z)) varies monotonically in the dust 
layer z/h^ < 2, and become almost constant in region outside the dust layer z/h^ > 2. 

In Figure 10 and Figure 11, the eigenfunctions of velocity are drawn in the case of 
a = 1cm and Pd/Pg = 1- Since the friction is sufficient to couple the gas and dust firmly, 
there is almost no difference of velocities between gas and dust. \vgz \ and |fgy| have maximum 
values at z/h^\ ~ 1.8, 1.3, respectively. As shown in Figure 15, the shear rate has a maximum 
value at z/h^ ~ 1.4. The velocity field of gas is drawn in Figure 14. Wc can sec the vortex 
clearly at z/hd = 1.5 which is near the point of the maximum shear rate. The unperturbed 
quantities interact with the perturbed quantities at the point of the maximum shear rate 
physically, and the perturbed kinetic energy of the y-component is supplied by the shear (see 
detailed discussion in Sekiya & Ishitsu (2000)). In Figure 12 and Figure 13, the eigenfunctions 
of velocities are drawn in the case of a = 100cm. The size of dust aggregates is too large for 
strong coupling, thus the velocity difference is large. As shown in Figure 11 and Figure 13, 
d\vgyi\/dz is discontinuous at z/h^ = 2, since there is no boundary condition that |fgyi| is 
differentiable. On the other hand, in Figure 10, |fgzi| is differentiable owing to two boundary 
conditions at z = 2/1^: Vg^i/i^^ — kU) and the pressure perturbation are continuous. 

The growth rate of the mode with the most unstable wave number as functions of the 
dust to gas density ratio on the midplane and the size of dust aggregate is drawn in Figure 
16. The growth rate increases as Pd/Pg increases. If pd is large, the thickness of the dust layer 
hd — Sd/pd is thin, and the velocity difference is large, thus the shear becomes large. Hence 
the dust layer becomes more unstable as the dust setthng proceeds. When Pd/Pg > 1, under 
the condition that Pd/Pg is constant, the growth rate has the local minimum at a ~ 1cm. 
When a > 100cm, the growth rate decreases significantly. The reason of the reduction of the 
growth rate at a ~ 1cm is caused by the friction discussed in Section 2. However the reduced 
growth rate is much larger than the Kepler frequency, and the shear instability is induced, 
dust aggregates are stirred up by turbulence. The reduction in the case of a > 100cm is 
not caused by the friction. The friction is so small that the unperturbed velocity is almost 
constant. Thus the shear is not so large to induce the Kelvin-Helmholtz instability. 

We compare the numerical result and the analytic result discussed in Section 2. We can 
expect that the characteristic length scale L is the thickness of dust layer. On the other hand, 
the analysis in Section 2 has no characteristic length scale. If we introduce the characteristic 
length scale into the analysis in Section 2, we may compare the result in this Section with 
the result in Section 2. Thus we consider k as 27r/L. When the stabihzation due to the 
inertia of dust is not effective, the maximum growth rate may be expressed as 27rC/ / L where 
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U is the characteristic velocity difference. 

The velocity difference between the midplane and the gas dominant region is expressed 
as the equation (57). In the small friction limit (^(pg + pa) ^ ^fe), we have 

5^;g~7^d(l + 7^d)^Vfc■ (80) 
In the strong friction limit, we have 

/Cd + i 

The thickness of the dust layer is 

Zd ~ Ed/pd- (82) 

With (28), (81), and (82), we obtain the condition of stabihzation: 

20 

—cm < a < 10cm, (83) 

where we assume T^d ^ 1 and fi = 0.5. As shown in Figure 16, the region where the growth 
rate is reduced owing to the friction corresponds to the condition (83), thus the numerical 
result is consistent with the condition discussed in Section 2. 

In the case of a ^ 1cm, the friction is very weak, and the reduction of growth rate due 
to friction is not effective. The friction coefficient A is proportional to in the Stokes 
regime. So we can estimate the growth rate /i: 

poca-^7^2(l + 7^d) (84) 



Thus the maximum growth rate is inversely proportional to in the case of a 3> 1cm, and 
proportional to TZ^ in the case of TZd < 1, or proportional to T^-d in the case of TZd > 1. 
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Fig. 6. — The distribution of the unperturbed dust density pdo i^^ft panel) and the unper- 
turbed velocity {right panel), where Pd/Pg = 1, a ~ l/im, 30cm, 100cm, — 0.5zd ■ We use 
a frame of reference moving with Kepler velocity. 
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Fig. 7. — The growth rate of the Kelvin-Helmholtz instabihty as a function of wave number 
in the case for a — 0.01cm, Pd/Pg = 0.1, where = The sohd hne shows the growth 

rate for the even mode. The dotted hne shows the growth rate for the odd mode. 
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Fig. 8. — Same as Fig.7, but in logarithmic scales. The solid hne shows the numerical result 
for the odd mode. The dotted hne and the dashed hne show the analytic solution (44) and 
(45) respectively. 




Fig. 9. — The amplitude {left panel) and the phase {right panel) of the eigenfunction Pi{z) 
in the case of a = 1cm and Pd/Pg = 1- 
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Fig. 10. — The amplitude {left panel) and the phase {right panel) of eigenfunctions Vgzi{z) 
(sohd curve) and Vd.zi{z) (dashed curve) in the case for a — 1cm and Pd/Pg = 1- The dashed 
curves are not seen because they overlap on the solid curves. 
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Fig. 11. — Same as Fig. 10, but for the eigenfunctions Vgyi{z) (solid curve) and Vdyi{z) (dashed 
curve). The dashed curves are not seen because they overlap on the solid curves. 




Fig. 12. — Same as Fig.lO, but for a = 100cm and Pd/Pg = 1- 
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Fig. 14. — The velocity field of gas in the case where a — 1cm and Pd/Pg — 1- 
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Fig. 16. — The contour of the peak growth rate of the Kelvin-Helmholtz instabihty as 
function of the dust to gas density ratio on the midplane and the size of dust. The growth 
rate is normahzed by the Keplerian frequency. 
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6. Discussion 

When the size of the dust grain a is smaller than 1cm, the stabilization by the friction 
is not effective. Thus the dust layer is expected to become turbulent, if the settling onto 
the midplane proceeds. When a > 100cm, however, the growth rate of Kelvin-Helmholtz 
instability is small, the settling may proceed without stirring up. In this section we discuss 
a possible path to the planetesimal formation through the gravitational instability. 

We compare the various timescales for dust grains in the disk as functions of the size 
of dust and the degree of the dust setthng onto the midplane. For simplicity, we assume 
that all dust aggregates have the same size. We consider the timescales of the dust setthng, 
Kelvin-Helmholtz instability, the dust growth, and the migration into the central star. 

The dust growth time can be calculated by 

_ 4 Pmat a 

f-grow — 7 T~! V^^J 

Js Pd dv 

where fs is sticking probability, pmat matter density, and Sv the averaged velocity dispersion 
of the dust aggregates. Sv is given by 

5v — max{5vi, SV2), (86) 

where 

Svi = ^ f^Vgt, (87) 

mn hydrogen mass, dust mass, z — z^, Vgt the thermal velocity of gas. 

In the laminar flow, we use the thermal velocity 5v — 5vi given by the energy equipar- 
tition of gas and dust. This is expected if the dust aggregate is small. If the dust aggregates 
are large and setthng, the dust aggregates collide with the other small dust aggregates which 
are not settling. In this case the relative velocity 5v might be determined by the settling 
velocity Sv2- 

The sticking probability is nearly unity in the case of the collision with small relative 
velocity. In the case of large relative velocity, the sticking probability decreases to zero 
(Poppe et al. 2000; Sekiya & Takeda 2003). However we use the sticking probability = 1 for 
simplicity, because the smaller sticking probability should not change the following diagram. 

In addition we use the equation (88) for large dust aggregates. In the laminar flow the 
flnal sizes of dust aggregates are 20cm at lAU (Nakagawa et al. 1986). In the case of large 



dust aggregates, we can not use the equation (88) to calculate the time scale of the growth of 
dust aggregates. We can use the equation (88) for small dust aggregates. On the other hand, 
if the dust layer is turbulent the equation (88) cannot be used. However we are interested 
in the boundary line between turbulent dominant state and settling dominant state. This 
boundary state does not depend on the growth time if the growth time of the dust aggregates 
is longer than the settling time. The settling for large dust aggregates is very fast, thus we 
can expect that the setthng time is shorter than the time scale of growth in the case of large 
dust aggregates. Thus we use the equation (88) for simplicity. The boundary hue should 
not change even if we use more realistic growth time. 

If a dust aggregates is small (a < 20cm), the settling time onto the midplanc is deter- 
mined by 

^settle = (^^) 

However if a dust aggregates is large (a > 20cm), it revolves around the central star decreas- 
ing the inclination gradually (Nakagawa et al. 1986). The decreasing timescale of inclination 
is order of stopping time l/Apg. Thus in the case of a > 20cm, the effective settling time is 
considered as decreasing timescale: 

^settle = iMPg- (90) 



The timescale of Kelvin- Helmholtz instability is 



Im[a;inaxJ 



The result is summarized in Figure 17. In the region of a < 0.1cm and TZd < 0.01, the 
growth of dust aggregates is the fastest process, thus the dust can grow up to a ~ 0.1cm. 
In the region of a < 0.1cm and TZ^ > 0.01, the Kelvin-Helmholtz instability is the fastest 
process. In the region of a > 0.1cm, the dust settling is faster than the dust growth, thus 
we should compare two processes: the settling and the Kelvin-Helmholtz instability. In the 
region of a 0.1cm, the settling velocity is very fast and the Kelvin-Helmholtz instabihty 
is not effective, thus the settling region expands into higher Pd/Pg region as the size of dust 
a increases. 

We assume that the turbulent flow become laminar if the timescale of Kelvin-Helmholtz 
instability is larger than those of other processes. If the density and the size of dust aggregates 
are Pd/Pg < 0.01, a — 10~^cm, the dust aggregates grow up to a ~ 0.1cm without the dust 
settling. When a > 0.1cm, the dust aggregates start setthng and the Kelvin-Helmholtz 
instability occurs, which develops turbulence in the dust layer. The dust aggregates may 
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grow owing to the random velocity field in the turbulent flow, thus the timescale of the 
Kelvin- Helmholtz instability will increase with the growth of dust. When the size of dust 
grain becomes sufficiently large, txH > Nettie; the dust aggregates begin to settle. Hence 
the dust layer may evolve along the boundary line between the turbulence region and the 
settling region in Figure 17. The density of dust in the disk may reach lOOpg at a ~ 10m. 

Pd ~ lOOpg corresponds to the order of the critical density for the gravitational insta- 
bility. However this argument is not conclusive, because the critical density may be larger 
than lOOpg. A more realistic critical density depends on a unperturbed dust density pro- 
file. Actually Yamoto & Sckiya (2004) conclude that the critical density for the constant 
Richardson number dust density distribution is 760pg and that for the Gaussian dust density 
distribution is 340pg. Thus the realistic condition of the gravitational instability seems to 
be depend on the circumstances. 

Ishitsu & Sekiya (2003) showed that the shear instability which grows slower than the 
Kepler timescale without the Coriolis and tidal forces is suppressed if these forces are taken 
into account. Their analysis is based on the strong coupling approximation but our analysis 
is not based on it, thus we cannot use their result directly. To take into account the Coriolis 
and tidal forces, we assume that the region txH > 27r/Qk should correspond to a stable 
region, where tKH is the timescale of the shear instability without the Coriolis and tidal 
forces, 27r/r2k is the Kepler timescale. We plot the stable region in Figure 18. This effect 
causes little change on a possible path. 

The growth timescale for meter-size particles is longer than the timescale of the migra- 
tion into central star because of the large radial vclocity(Adachi, Hayashi, & Nakazawa 1976; 
Weidenschilling 1977). This is correct when the dust density is smaller than the gas density 
(T^d < !)• The radial velocity in the case of TZ^ > 1 is slower than the radial velocity in the 
case of T^d < 1 (Nakagawa et al. 1986). The radial velocity of dust is expressed as 

2r 

Thus the migration timescale tmigration = ''"/I'i^drl is 

_ (l + (l+7^d)^^> .gg. 

In the limit of TZ^ ^ l/F, the migration timescale is proportional to (1 + T^d)^- In the case 
where TZ^ — 0.01 and a — 100cm, the migration timescale is about 10^ year at 1 AU. In 
contrast, if TZ^ — 100 and a — 100cm, the migration timescale is about 10^ year. Thus the 
dust aggregates can avoid migration onto their central star and lead to the gravitational in- 
stability, if the growth timescale in the turbulent flow is shorter than the migration timescale 
10^ year. 
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Fig. 17. — The dominant processes in the disk of various parameters. The Kelvin- Helmholtz 
instabihty region is denoted by the crosses. The thin arrows mean the variation of a and 
Pa/ in the region where the sedimentation or the growth is dominant process. The block 
arrows mean the growth in the turbulent flow. The bold arrow is the possible path to the 
planetesimal formation through the gravitational instability. 
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Fig. 18. — Same as Figure 17 but we take into account the stabilizing effect of the Coriolis 
and tidal forces. Filled circles correspond to the models that are expected to be stabilized 
by Coriolis and tidal forces. 
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7. Conclusion 

In this paper, we have investigated the Kelvin-Helmhohz instabihty for the dusty fluid 
taking into account of the effect of relative motion between gas and dust. 

To clarify the physics of Kelvin-Helmholtz instability in the dusty fluid, we have analyzed 
the linear stability of the flow whose velocity profile is the simple Heaviside step function. 
We have found the stabilizing effect due to friction. In the case of pd/Pg ^ 1, the reduced 
growth rate is 

^(Pd/Pg)-'/% (94) 



Im[a;s 

where Im[a;sc] is the growth rate of the Kelvin-Helmholtz instability for simple fluid without 
dust aggregates. 

We have considered the case of dusty layer between two gaseous layers. Because of the 
symmetry, there are the even mode {Pi{z) = Pi{—z)) and the odd mode {Pi{z) = —Pi{—z)). 
We have found that the growth rate for odd modes is more unstable than that for even modes. 

We have analyzed the linear stability in a more realistic protoplanetary disk where the 
dust density distribution is a sinusoidal density distribution. In the case of a < 0.1cm, we 
found the stabilization due to the friction. This result is consistent with the result in the 
case of the step function velocity profile. However the timescale of the stabihzed Kelvin- 
Helmholtz instability is shorter than that of settling, thus the dust layer becomes turbulent 
before the gravitational instability sets in. In the case of a ^ 1cm, since the shear is small, 
the Kelvin-Helmholtz instability is stabilized sufficiently. In this case, if the dust aggregates 
begin to settles, it is expected that the dust settles without stirring up and may reach the 
critical density. 

We have compared the timescales of the dust growth, the dust settling, and Kelvin- 
Helmholtz instability, and proposed a scenario in which the dust layer proceeds along the 
boundary between the turbulent region and settling region in Figure 17, and the dust layer 
reaches the order of critical density of gravitational instability at a ~ 10m. In this way, the 
planetesimals may be formed through gravitational instability. 

To justify this scenario, we have to do the detailed numerical calculation of the growth of 
dust aggregates in the turbulent flow. In this paper, we neglected the size-distribution of dust 
grains whose effect should be studied. We expect the density of dust layer becomes critical 
density for gravitational instability with the help of turbulence. However the dust grains in 
the turbulent flow have the velocity dispersion. Thus the critical density for gravitational 
instability possibly changes owing to the increased velocity dispersion of dust grains. We 
have to study the criterion of gravitational instability in more detail when the dust grains 
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have the velocity dispersion in turbulent flows (c.f. Weidenschilling 1995). In subsequent 
papers, we plan to study these effects. 

The authors thank the anonymous referee for valuable comments and suggestion that 
improved the paper. This work is supported by the Grant-in- Aid (No. 157401 18, 16077202, 
16244120) from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) 
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A. The effect of the velocity dispersion of the dust grains 

We consider the velocity dispersion of the dust grains. For the sake of simphcity, we 
assume that dust component obeys the equation of state of an ideal gas, 

^'d = 4pd, (Al) 

where Pd is the dust pressure, Cd the dust isothermal sound velocity. The equations governing 
this system are 

V • vg = 0, (A2) 
^ + V(pdVd) = 0, (A3) 



Pg + (vg • V)vg^ = -VP, + ^PgPd(vd - Vg), 



(A4) 



Pd + (vd • V)vd^ = -CdVpd - MPd(vd - Vg). (A5) 

We consider a Heviside step function velocity profile, and uniform density as the unper- 
turbed state: 

[/_ for z < 



We restrict ourselves to the mode in the direction of the y-axis. The hydrodynamic equations 
are given by 

ikv,yi{z) + ^{z)^0, (A7) 
az 

-i{uj - kU±)pdi{z) + ikpdoVdyi{z) + Pdo-^(^) = 0> (A8) 

-ApgpPdo{vdxi{z) - Vgxi{z)) - ipgpVgxi{z){uj - kU±) = 0, (A9) 

-Apgopdo{vdyi{z) - Vgyi{z)) - ipgoVgyi{z){uj - kU±) + ikPi{z) = 0, (AlO) 

dP 

-Ap^Pdo{vdzi{z) - Vgzi{z)) - ip^Vgzi{z){uj - kU±) + -^{z) = 0, (All) 

-^PgoPdo(^gxi(2;) - Vdxi{z)) - ipdQVdxi{z){'^ - kU±) = 0, (A12) 

-Ap^Pdo{vgyi{z) - Vdyi{z)) - ipdoVdyi{z){uj - kU±) + ikc\pd\{z) = 0, (A13) 

-Ap^QPdo{v^zi{z) - Vdzi{,z)) - ipdoVdzi{,z){uj - kU±) + ^d-^l^) ^ 0- (A14) 

We transform the physical quantities to nondimensional ones by taking kr^. p^^kr^ , and 
{kU)~^ as units of the length, the mass and the time, respectively, where U = ^7+ — 
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We rearrange the equations to obtain the matrix form of the hnear equations: 

n_ 





( Pii^) \ 
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dz 
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^d drag 
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Vdzl{z) 

(A15) 


/ 





where 



N = -cliTZd-itd, 



rag 



Uo)) + (1 + Ted - iidrag(^ - t/o))(u; - Uo) 



We define the vector and the matrix as v and A respectively. 



(A16) 
(A17) 



V is expressed as the hnear combination of exp(Ai2)vj where i — 1,2,3,4. Aj and Vj are 
given by 



(Ai, A2, A3, A4) — (1, —1, Aa, — Aa), 

(Vi,V2,V3,V4) = 

// 1 \ / 

(l-ttdrag(^-t/o)) 



D 



D 



1 \ 

(l-itdrag(t^-i^o)) 



D 







D 



( 



1 



i{A2-l)t, 



drag 



■^a^drag 

Ted 



1 



i(Ag-l)tdr, 
Aa^drag 



(A18) 



/ 



where 



Aa = 



Cd'x/^drag 



(A19) 



(A20) 



We cannot determine the sign of Re(Aa) if we do not know lo. However if we assume 
Re(Aa) > to calculate uj, we can obtain the consistent solutions. Thus we assume Re(Aa) > 
0. (In addition, in the case of Cd = and tdiag = 0, wc can show Aa > 0.) 

We seek the solutions which does not diverge in the limit of \z\ 00. In the case where 
z > ^, the coefficients of Vi and V3 are zero. In the case where z the coefficients of V2 
and V4 are zero. 

Thus we have the conditions: 

+ 7ed) + (-| + c^)idrag)^^gz+ +u;)pd+ 



rag 



^drag( — 1 + Aa+) 



0, (A21) 
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p (|+a;)(l + 7ed-i(|+u;)idrag)^^gz- i(| + a^)pd- „ 

-Tl- i \ -j r -. y , — U, [AZO) 

I + {2+ ^)Wg Wgl-l + Aa_j 

^d^d.- 2i7^d^;g.- i(|+a;)pd- ^ .^24) 

Aa- (2z+(l + 2u;)Wg)Aa- -l + A^ ' ^ ' 

The boundary conditions are determined by the continuity of pressure and Vgzi/i^^ — kU), 
we have 

Pi+ = Pi-: (A25) 

Pdi+ = pdi-, (A26) 

^ = (A27) 



^+2 

These equations have a nontrivial solution, if the determinant of the matrix of the hnear 
equations is zero, thus we have the dispersion relation: 

(le^drag'^^ + 16i(2 + 7ed)tdragC^' " 16(1 + 11^)^^ - 4^(7^d - 2)tdragO^ - 47ed - 4 - ^drag') 

x((2u; - l)(tdrag(2cU - 1) + 2i)Aa- - {2u + l)(idrag(2a; + 1) + 2i)Aa+) 

+8nd (4a;2 - l) Cd^drag = 0. (A29) 

In the limit of Cd 0, this dispersion relation corresponds to the dispersion relation in 
Section 2. Thus we can expect that the analysis in Section 2 is valid although Vdzi{z) / (u—kU) 
has the discontinuity at the boundary. Vdzi{z)/{u! — kU) becomes discontinuous at the 
boundary in the case where Cd = although that is continuous for any case where Cd 7^ 0. 

To understand the reason of this transition from continuity to discontinuity in the limit 
of Cd — > 0, we investigate the velocity eigenfunction of dust. The velocity perturbation of 
dust is given by 

vazi{z) = Ci exp(-2;) + C2 exp(-Aa2;), (A30) 

where Ci and C2 are constant. This is shown in Figure 19. The dotted line corresponds to the 
case where Cd — > 0. The correction to the Cd — > case is C2 exp(— Aa-z). i'd2i(-z)/(<^ — kU) is 
continuous because of this correction, that is, if we consider only Ci exp(— 2;) , v^zi i^) / (uj—kU) 
is discontinuous at the boundary. This correction term becomes significant in the region 
where \z\ < A. Figure 20 shows that A diverges in the limit of Cd 0, thus the width of the 
region for the correction decreases to 0. Hence the correction term disappears effectively. In 
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Fig. 19. — The amplitude of the eigenfunction Vdzi{z) in the case where Cd — 1 (sohd hne) 
and the amphtude of the eigenfunction v^ziiz) — Ciexp(— in the case where Cd — 
(dotted hne). 



this way we can understand the transition from continuity to discontinuity in the hmit of 

In contrast, Figure 21 shows that the growth rate converges to that of ca = in the 
hmit of Cd — >^ 0. Thus our resuh from the analysis with discontinuous velocity perturbation 
for dust with Cd — can be regarded as a rigorous result for the case where Cd is very small. 

In the limit of Cd — > oo, the growth rate approaches to that of simple Kelvin-Helmholtz 
instability. As shown in Figure 22 and Figure 23, the stabilization by the friction lessens as 
Cd increases. 
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Fig. 21. — The the growth rate as a function of c^. The sohd Une shows the growth rate 
where — 1 and the dotted hne shows the growth rate where ca = 0. 
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Fig. 22. — The contour of growth rate as a function of idrag and TZ^ in the case where ca = 1. 
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Fig. 23.— Same as Fig.22 but for Cd = 0.1. 
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B. The most unstable mode 

We show that the solution of the equation (27) which has the largest real part is real. 

We define the left hand side of the equation (27) as /(cD). We can show /(O) < 0, 
/(oo) = oo, /(— oo) = oo. It follows that the equation (27) must have at least one positive 
solution Xi and one negative solution X2- If two other solutions are real, it is true that the 
solution that has the largest real part is real. In following discussion, we consider the case 
that the equation (27) has two complex solutions. 

The complex solutions are 2:3 and its complex conjugate because all coefficients are 
real numbers. We can rewrite the equation: 

(o;^ - 2luj - n?) {uJ^ - + n^) = 0, (Bl) 

where 

m = \/—XiX2i (B2) 

n = |x3|, (B3) 

1^x1+ X2, (B4) 

7 = Re(x3), (B5) 
This equation is equivalent to the equation (27), we obtain the equations: 

-(^ + 7) = ^^, (B6) 



rag 



-m' + n' + 4/7 = —2 ^ , (B7) 

Mrag 

7 2 I 2 "I" ^d^drag fDo\ 

—In +mj — -, (B8) 

^drag 



drag 

Prom the equation (B6) and (B8), we give 



7 -2 + 7^d + m^(2 + 7^d) 



\-_ 

[m? + n2)tdrag 



,= Z1±]1±ZJ1^1±M. (Bll) 
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In the case of T^-d < 2, from the equation (Bll), we can see 7 < . 

In the case of TZ^ > 2, we give the inequality / < from the equation (BIO). Here after 
we use the reduction to absurdity to prove 7 < 0. We suppose 7 > 0. Prom the equation 
(B7), it can be seen that — > 0. Prom the equation (B9) we can see rn^n^ > 1. Then 
we can see > 1. On the other hand, from the equation (Bll) we give the inequahty 



7^d + 2 



(B12) 



It can be seen that 



7^d + 2 



< 1. This imphes that < 1, a contradiction. Thus we obtain 



the inequahty 7 < 0. 



This means that the real parts of the complex solutions are always negative. Thus the 
solution that has the largest real part is not complex but real: the most unstable mode is 
the exponentially growing solution. 



